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Abstract 

What basic processes generate hierarchy in a collective? The Bonabeau model pro- 
vides us a simple mechanism based on randomness which develops self-organization 
through both winner /looser effects and relaxation process. A phase transition be- 
tween egalitarian and hierarchic states has been found both analytically and numer- 
ically in previous works. In this paper we present a different approach: by means of 
a discrete scheme we develop a mean field approximation that not only reproduces 
the phase transition but also allows us to characterize the complexity of hierarchic 
phase. In the same philosophy, we study a new version of the Bonabeau model, de- 
veloped by Stauffer et al. Several previous works described numerically the presence 
of a similar phase transition in this later version. We find surprising results in this 
model that can be interpreted properly as the non-existence of phase transition in 
this version of Bonabeau model, but a changing in fixed point structure. 
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1 Introduction 



It is usual, in sociological works, to describe how global behavior appears, in 
many levels of social activities [1]. Before that, it is more fundamental to un- 
derstand in which way citizens gather [2] : since in a little collective every single 
seems to play the same status, in big societies diversity appears [3] . Hierarchic 
dominance and hierarchic stratification has been studied with several different 
approaches [2,4,5]. As long as these matters can be considered as many-body 
dynamical systems, they have attracted the attention of physicists in the latest 
years. The emergent area of Sociophysics involves those social complex sys- 
tems, dealing with many different social situations with a statistical physics 
approach [8,1,3]. 

In this terms, a simple and fruitful model of diversity generation has been 
proposed by Bonabeau et al. [7,6]. Related to this model, some modified ver- 
sions have been proposed, as the Stauffer et al. [10] version, or the one from 
Ben-Naim and Redner [16] (this later one has been solved analytically). 
The purpose of this paper is double: first, by means of a discrete mean 
field approximation, we reproduce analytically the numerical results found 
by Bonabeau. We discover a non trivial complex structure in the hierarchy 
generation path. After this, we apply the same scheme to Stauffer version 
[10,9], widely used as a model of hierarchy generation [10,9,13,14,11,15], in 
order to obtain analytical evidences of its numerical behavior. 



2 The Bonabeau model 

The Bonabeau model [7,6] has been proposed as a simple model showing 
self-organization to explain hierarchic dominance in Ethology. With subtle 
modifications it has been reallocated in Sociophysics area as a model of social 
stratification [9,10,11,16]. It starts from a community composed of N agents, 
randomly distributed over a regular lattice L x L, that is to say, with a pop- 
ulation density p = Each agent i — 1,2, N is characterized by a time 
dependent variable hi(t), the agent fitness, that from now on we will name 
status. Initially all agents share the same status hi(t = 0) = 0: the so-called 
egalitarian situation. System dynamics are: 

(1) Competition with feedback: an agent % chosen randomly moves in a four 
nearest neighbor regular lattice (Newmann neighborhood). If the target site 
is empty, the agent takes the place. If it is already occupied by an agent j, a 
fight occurs. The attacker agent % will defeat agent j with some probability: 
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1+exp {ri(hj(t) - hi(t))) ' 
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Where 77 > is a free parameter. If i wins, he exchanges positions with j. 
Otherwise, positions are maintained. After each combat, status hi(t) are up- 
dated, increasing by 1 the winner's status and decreasing by F the looser's 
status. Note that F is a parameter of the system that weights the defeats such 
that F > 1. The case F — 1 will be the symmetric case from now on. In the 
asymmetric case, F > 1, the fact of loosing will be more significative for the 
individual status than the fact of winning [9]. 

(2) Relaxation: a natural (Monte Carlo) time step is defined as N movement 
processes (with or without combat). After each time step all agents update 
their status by a relaxation factor (1 — //), such that < /i < 1; this effect is 
interpreted as a fading memory of agents. 

Notice that competition rule (1) is a feedback mechanism: status differences 
hj(t) —hi(t) drive the future winning/loosing probabilities of agents i and j. If 
agent % wins/looses it's winning/loosing probability increases afterwards. This 
mechanism amplifies agent inhomogeneity. On the other hand, relaxation rule 
(2) drives the agent status hi(t) to equalize: status differences are absorbed 
and toned town. The balance between both mechanisms generates asymptotic 
stability on hi(t). Common sense would lead us to expect low fights when the 
agent's density is low, so that the relaxation mechanism would overcome and 
egalitarian situation would prevail (hi = hj = 0,Vi, j). But if the system pos- 
sesses high agent's density, the rate of fights would increase, and competition 
mechanism will prevail, leading the system to non neglecting inhomogeneities. 
The balance of these two mechanisms is crucial at a given density p. Simu- 
lations ran by Bonabeau et al.[7,6] show how this compromise between both 
effect bring about a phase transition at a critical density, between egalitarian 
societies for low densities and hierarchical societies for high densities. 

A natural measure for the status diversification is the standard deviation of 
its stationary distribution {/i*}i=i,...,jv- However, in [9] another measure is pro- 
posed: this is the standard deviation of stationary probability distribution 
{Pij}i=i,...,N, defined as: 



This choice turns out to be more suitable , as long as it is a bounded parame- 
ter: a G [0, 1]. It works as an order parameter of the system: for densities lower 
than the critical, every hi are equal and therefore all P^ too, then o = 0. For 
densities bigger than the critical, status are different and therefore probabili- 
ties are also different: this leads to a non zero value of a. 




(2) 



3 



3 Mean field approximation in the Bonabeau model 



In order to tackle the system in a mathematical way, we will obviate spatial 
correlations, reinterpreting p as the probability of two agents combat. In the 
spatial correlated model this is equivalent to a random mixing of the agent's 
positions every time step. Therefore, at each time step, an agent % will possess: 

(1) Probability 1 — p of no combat. In that case, the agent will only suffer 
relaxation. 

(2) Probability p of leading a combat (with probability l/(N — 1) the attacked 
agent will be j). In that case: agent % will increase, in average, its status by 
one with probability Pij(t), and will decrease by F its status with probability 
1 — Pij(t). Relaxation will also be applied in this case. 

The model is then described by an N equation system {i = 1,...,N) of the 
following shape: 

hi(t + 1) = (1 - p)(l - AW) + f; {Py(f)(fc(*) + 1) + (1 - Pij(t))(hi(t) ~ F)}-(3) 

In order to analyze the system let's start with the simplest version N = 2. 
Having in mind that for two agents Pi 2 (t) = 1 — ^2i(0> the system (3) will 
reduce to: 

h(t + 1) = (1 - A*)/n(t) + p(l - p){Pi2{t){\ + F)-F} 

h 2 (t + 1) = (1 - p)h 2 (t) + p(l - p){\ - P 12 (t)(l + F)}. (4) 



In the egalitarian phase (below the critical density), the fixed point (h{, h^) of 
the two-automata system 4 will have stationary status of the same value, say 
h\ = h* 2 . In order to find the fixed point of the system we can define the mean 
status of the system as: (h(t)) = (hi(t) + h 2 (t))/2. The system turns into: 

(h(t + 1)) = (1 - p)(h(t)) + P(1 "^ (1 " F) , 
with a fixed point: 

</>•> = ^-f- F \ (5, 



always stable (this can be interpreted as the system's energy, which is con- 
served) . 



4 



Notice that, in the egalitarian phase (below critical density), we have h\ = 
h\ = (h*). In order to check out stability of ({h*), (h*)) we compute the 
Jacobian matrix of the system, evaluated in that fixed point: 



/ 



I - A A 



\ 



J = (l-li) 



(6) 



V 



A I- A 



) 



where: 



A = 



(7) 
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By stability analysis we conclude that egalitarian phase is stable as long as: 



From numerical iteration of the equations for two-automata, we present in 
figure 1 control parameter p versus order parameter a, at stationary situation, 
for different values of parameters F, 77 and p. Critical density values p c agree 
with (8). For the symmetric case (F = 1) with 77 = 1 and p, — 0.1 (squares) 
we obtain a critical density p c ss 0.11. We'll take this particular case as the 
reference case from now on. 

In order to understand the dynamics of the fixed points of the two-automata 
system, we apply the following change of variables: h* = h\ — h\. The fixed 
points of the system become now the solutions h* of: 



As to the case of reference, notice that for each value of p the system, as a 
fact of symmetry, has two fixed points: (h\, —h\), (h^, —hi). In figure (2) we 
represent a case below transition (p < p c ), with a singular solution, and a 
case above transition (p > p c ) with three solutions: this results in the solu- 
tion bifurcation of h\ and h* 2 in figure 3. In the egalitarian phase we have: 
h\ — h* 2 — (h*) = 0. At p = p c a bifurcation occurs, from which we have 



In the successive figures we can observe what effects produce a variation of 
parameters F, 77 and p on status stratification, as for the reference case (fig- 
ure 3): (1) Increasing the asymmetry F (figure 4) decreases p c and grows up 
inequality. (2) Increasing the relaxation p (figure 5) increases p c and diminish 
inequality. (3) A decrease of 77 (figure 6) increases rho c and has no effect on 



P< Pc 



2p 



(8) 



77(1 -p)(l + FY 




h{ = -h\ ± 0. 
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inequality. 



We now apply the same philosophy to the general system 3. We may define the 
mean status of the system as: (h(t)) = jj X^Li hi(t). Now due to the fact that 
Pij + Pji = 1 and then Y^iLi X^lij^i P%j — N ^ N ~^> , we can sum and normalize 
the N equations, obtaining: 



{h(t+1)) = (l -, ){m - ^-^ F -») , 

whose fixed point is independent of the number of automata and agrees with 
the very first result given at (5) in the case of two-automata (the system's 
energy). Again we get that (h{, h^, ■ ■ ■ , h* N ) with h\ — h\ — . . . — h* N — (h*) 
is a fixed point of the system whose stability determines the transition. The 
Jacobian matrix of the linearized system, evaluated at the fixed point, is: 



J 



I- A 



N-l 



\ N-l 



N-l 

I- A 



N-l 



I- A 



a circulating matrix [12] where A = anc i its eigenvalues being: A = 

(1— A) — with multiplicity TV— 1, and A = (1— //) with multiplicity 

1. The egalitarian phase is therefore determined by: 



P < Pc 



77(1 — n)N(l + F) 



This result is on agreement with the particular case of two-automata (N — 2), 
and if iV >> 1 we have: 



Pc 



4/x 



r/(l-/i)(l + F) 



(9) 



Notice that as long as < p c < 1, the phase transition will only occur under: 
77(1 + F) , . 

"<irwTFy (10) 



In figure (7) we represent, for iV >> 1, the parameter space, where we distin- 
guish the zone where the egalitarian-hierarchical transition is allowed. 
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4 Additive relaxation 



Bonabeau et al. in their seminal paper [6], proposed an additive relaxation as 
an alternative to the multiplicative relaxation developed above. That additive 
relaxation updates the status as it follows: 

hi — > hi — /itanh(/ij). (11) 



They developed then a mean field approximation, based on stochastic differen- 
tial equations, where they found the egalitarian-hierarchical phase transition. 
We can of course apply, in the discrete model that we propose, this additive 
mechanism of relaxation. In the two-automata system, mean field equations 
reduce to: 

h(t + 1) = h^t) + p(P 12 (t)(l + F)-F) - /itanh(/n) 
h 2 (t + 1) = h 2 (t) + p(l - P 12 (t)(l + F)) - pttmh(h 2 ). 

Following our previous steps, it's quite easy to deduce that the fixed point 
hi = h* 2 is stable as long as p < p c = V ^+ F ) ■ This result is on agreement with 
those of the continuum model proposed by Bonabeau et al. [6]. 

A simplified model has been recently proposed [16] by Ben-Nairn and Redner, 
in order to obtain analytical evidences of the phase transition which is observed 
in the Bonabeau model. In their version (which is highly inspired on the 
Bonabeau model) the relaxation process is additive (though it is not a function 
of the status but simply a fixed constant hi — > hi — 1), and the competition 
mechanism is not stochastic but deterministic (except for the situation where 
both agents have the same status). A phase transition between both regimes 
(equality/hierarchy) is then found analytically. 

Comparing additive to multiplicative relaxation in the Bonabeau model, we 
must say that additive works worse than multiplicative: status differences 
grow excessively and then, computing limitations are exceeded (due to the 
exponential explosion). 

Anyway, if we introduce a new parameter Q > (instead of increasing by 
one the winner status, we increase it by Q, so that we can tune both winning 
and losing effects Q,F), the system doesn't explode numerically with a good 
tuning of the parameters. In figure 8 with N = 10, we took F = 0.7, Q = 0.7, 
p = 0.0001 and rj = 0.001. As we see, status reach values of five order of 
magnitude, what we think is not realistic. 

The hierarchy scheme in the multiplicative relaxation model (figures 9-13) 
develops much more complexity than the additive one. An increase of density, 
above critical one, stills generates hierarchy, as common sense would have 



7 



dictated us. This fact is traduced by a periodic fixed point coordinates h* 
splitting, even after the phase transition. This hierarchy growing is not trivial; 
something that has for sure been unnoticed for the moment. In the additive 
relaxation model the hierarchical structure is simple, it doesn't change with 
p at hierarchical phase. Instead of that, there is a fixed point coordinates h* 
splitting at p c , and dynamical evolution in hierarchical phase is poor. 



5 Mean field in the Stauffer version 



The model developped by Stauffer et al. comes from Bonabeau's. It was firstly 
introduced to carry out the supposed lack of transition of the previous model, 
discussed numerically in [10,9,11]. In this version, the free parameter 77 is now 
exchanged with the order parameter a, such that 



Pll{t) l + e W (a(t)(h j (t)-h i (t))Y (12) 



This modification somehow introduces a dynamical feedback to the system: 
probability of winning/loosing is directly related to the global inequality of 
the system, therefore, depending on the climate's aggressiveness of the system, 
agents will behave more or less aggressive themselves. 

Just as in the case of Bonabeau's model, by introducing a mean field we expect 
to find and reproduce analytically the phase transition that is proclaimed in 
the literature. 



In the case of two-automata system, we develop the mean field equations hav- 
ing in mind that: 



(1) Each automaton updates at each time step with the same dynamics that 
the Bonabeau model. The probability calculation is different though (77 1— > a). 

(2) At each time step, variable a is updated: the order parameter is now a 
dynamical parameter of the system and therefore evolves with it. 



We once again redefine h — h 2 — h\. With this change of variable, we pass 
from a three equation system (hi, h 2 , er) to a two equation system (h, a), with- 
out lack of generality, because hi and h 2 are related through the mean status 
(system's energy). The system equation is therefore: 
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h(t + 1) = (1 - n)h(t) + p(l + F) 

1 1 , 

a(t + l) 



1 + exp (a(t)h(t)). 



1 + exp (a (t)h(t)) 2 1 ' 

It is likely to expect the same qualitative results of the stability analysis of 
this system than those about Bonabeau's model, that is to say, loss of stability 
of the egalitarian regime (i.e. the fixed point h*(— h\ — hi) = 0, a* = 0) would 
become unstable at some critical density p c ). 



The Jacobian matrix of the linearized system is, in general: 



J 

where: 



1 - p + ABa* 2ABh* 
-a*A -h*A 



A = exp(a*h*)/(l + exp(a*/i*)) 2 , 
B = p{l-n){\ + F). 

Evaluating J in the egalitarian fixed point (h* = 0;cr* = 0), we find that 
eigenvalues of J are: 



\ 1 = < 1,A 2 = 1 -n < l,Vp (13) 

We find that egalitarian zone is stable for all densities. How come a phase 
transition can then occur, as is presented in many previous model simulations 
[10,9,13,14,11,3,15]? How come hierarchical situation can be achieved starting 
from equality, if the egalitarian zone is always stable? The key of the dilemma 
is set on the simulation methods that have been applied until now. Figure 
(14) shows stationary values of order parameter a vs. control parameter p. 
For each p, we take random initial conditions for status and a (these can be 
both zero or non-zero). The figure shows stationary value a = for all initial 
conditions, below a certain density. Above it, we have stationary values of a 
being zero in some cases (we remain in the egalitarian zone) and non-zero in 
others (hierarchical zone). 

Numerically we observe that the system has one stable fixed point below 
the critical density p c (h* = 0; a* = 0) , and two stable fixed points above 
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(h* = 0; a* = 0, h* ^ 0; a ^ 0). The egalitarian zone is therefore always stable 
(Vp). At p c , a saddle- node bifurcation occurs, and brings about the hierar- 
chic branch. Notice that the stability scheme is totally different from what 
we founded in the Bonabeau model: while in that model, equality-hierarchy 
transition was generated across Pitchfork bifurcation, due to loss of stability of 
the egalitarian regime, in the Stauffer version the egalitarian regime is always 
stable, but here at p c a saddle-node bifurcation takes place. The stable branch 
of this bifurcation is related to the hierarchic regime, and the unstable branch 
(not drawn in figure 14) plays the role of frontier between the two domains of 
attraction. 

Depending what initial conditions we give to (hi, a), (i.e. depending in which 
domain of attraction we start), the system, above p c , will evolve towards the 
egalitarian domain or the hierarchic one. 

In figure 15 we can understand how this stability is developed. Below p c the 
system has only one fixed point (triangles), indeed stable due to (13): every 
set of initial conditions h, a will evolve towards h* = 0, a* = (egalitarian 
zone). Above p c the system has three fixed points (circles), moreover, from 
(13) and figure 14 we know that upper and lower fixed point are stable and 
characterize both egalitarian and hierarchic zones, this leads to an unstable 
fixed point between them, performing the frontier. If the initial conditions be- 
long to the egalitarian domain, the system will evolve towards an asymptotic 
egalitarian state. On the contrary, if the initial conditions belong to the hier- 
archic domain, the system will evolve to an asymptotic hierarchical state. 

If in the model simulations, we set egalitarian initial conditions, the result 
would be the "absence of transition". But if we fix some other initial condi- 
tions, this could lead us to interpret the results as "existence of transition". In 
simulations made by Stauffer [3] they say "for the first ten Monte Carlo steps 
per site, a — 1 to allow a buildup of hierarchies" . This fact probably allows 
initial fluctuations develop so that initial conditions will be in the hierarchic 
domain. 



6 Conclusions 

The Bonabeau model has been criticized [10,9,11] in the last years. In this 
paper we revisit Bonabeau model in order to obtain analytical evidences that 
give clear proof of the phase transition that the system shows. We obtain, 
in the model with multiplicative relaxation, a high complex structure of the 
hierarchical regime, a fact that we think deserves an in-depth investigation. 
The Stauffer version, which is an alternative to Bonabeau model, proposed 
by Stauffer et al. [10,9,11], is the base of recent works, basically focused on 
simulations [13,14,15]. In this paper we tackle this version with the same phi- 
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losophy applied in the Bonabeau model. Surprisingly, this one doesn't show a 
phase transition in rigor, as far as there is no sudden growth of hierarchy if 
we start from equality. 
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p 



Fig. 1. Equality- hierarchy phase transition with control parameter p and order pa- 
rameter a, in the mean field model of two automata (N = 2). From left to right we 
represent the transition derived from iteration of the two automata equation system 
4, for different values (F; n;rj): (2;0.1;1) diamonds, (1;0.1;1) squares, (1; 0.1; 0.5) 
left triangles, y (1;0.3; 1) circles. The continuous lines are just guides for eye. 



Fig. 2. The two automata system switches from having one fixed point 
h* = h\ — h\ = for densities p < p c , to three fixed points when p > p c : 
h* = {0, +a, — a} (first one unstable and the rest stable), equivalent to the three 
fixed points of the two automata system, say {(0,0), (+a/2, —a/2), (—a/2, +a/2)}. 
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p 



Fig. 3. Bifurcation of stationary values h\ and h\ of the reference case (sym- 
metric case F = 1 with rj = 1 and p = 0.1). At the egalitarian zone we have: 
h* = h\ = (h*) = 0. At p = p c a bifurcation occurs, from where, due to symmetry, 




Fig. 4. Bifurcation of stationary values of h\ and h* 2 for values F = 2, p = 0.1 and 
77 = 1.0. 




Fig. 5. Bifurcation of stationary values of /i^ and h* 2 for values F = 1, p = 0.3 and 

T) = 1.0. 
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p 

Fig. 6. Bifurcation of stationary values of h\ and h\ for values F = 1, \i = 0.1 and 
r] = 0.5. 




3 4 

TV(UF) 



Fig. 7. Parameter space: delimitation, for N » 1, of the regions where transitions 
can be wether achieved or not achieved. 




P 



Fig. 8. Bifurcations of stationary values of fixed point components (system 3) de- 
pending on p, for the symmetric case of reference, in the additive relaxation model, 
when N = 10. 
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Fig. 9. Bifurcations of stationary values of fixed point components (system 3) de- 
pending on p, for the symmetric case of reference when N = 3. 




Fig. 10. Bifurcations of stationary values of fixed point components (system 3) 
depending on p, for the symmetric case of reference when N = 4. 




Fig. 11. Bifurcations of stationary values of fixed point components (system 3) 
depending on p, for the symmetric case of reference when N = 6. 
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p 

Fig. 12. Bifurcations of stationary values of fixed point components (system 3) 
depending on p, for the symmetric case of reference when N = 8. 




Fig. 13. Bifurcations of stationary values of fixed point components (system 3) 
depending on p, for the symmetric case of reference when N = 10. 




Fig. 14. a vs p in the 2 automata system. Initial conditions for each p are taken 
randomly. This leads to stationary values of a being zero (egalitarian zone) or non 
zero (hierarchical zone) depending of those chosen initial conditions. 
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Fig. 15. Crossover from one stable fixed point of the system at p < p c (triangles) to 
two stable fixed points -egalitarian/hierarchic- (circles) and an unstable fixed point 
as the delimiting branch between both domains of attraction. 
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